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ABSTRACT 

Measurements of solar wind turbulence reveal the ubiquity of discontinuities. In 
this study, we investigate how the discontinuities, especially rotational discontinuities 
(RDs), are formed in magnetohydrodynamic (MHD) turbulence. In a simulation of 
the decaying compressive three-dimensional (3-D) MHD turbulence with an imposed 
uniform background magnetic field, we detect RDs with sharp field rotations and little 
variations of magnetic field intensity as well as mass density. At the same time, in 
the de Hoffman-Teller (HT) frame, the plasma velocity is nearly in agreement with 
the Alfven speed, and is field-aligned on both sides of the discontinuity. We take one 
of the identified RDs to analyze in details its 3-D structure and temporal evolution. 

By checking the magnetic field and plasma parameters, we find that the identified 
RD evolves from the steepening of the Alfven wave with moderate amplitude, and that 
steepening is caused by the nonuniformity of the Alfven speed in the ambient turbulence. 


1. INTRODUCTION 


Since the beginning of the space age, the solar wind has been regarded as an excellent natural 
laboratory for studying the plasma turbulence. The endeavoring research over decades has revealed 


that the discontinuities are ubiquitous in the solar wind 

(Burlaga|l969, 1971, Smith 1973 

Solodyna 

et al. 1977, Tsurutani & Smith 1979; Neugebauer et al. 1984. Lepping & Behannon 1986 Tsurutani 

et al. 1994, 1996 Lee et al. 1996 Horbury et al. |2001 

Knetter et al. 2004, Tsurutani et al. 

2005] 

Neugebauer 2006 Vasquez et al. 2007; Li 2008; Lin et al. 2009. Sonnerup et al. 2010 

Teh et al. 

2011. Haaland et al. 2012 Malaspina &; Gosling 2012 

Wang et al. 2013; Paschmann et al. 

2013). 


These discontinuities appear as large and rapid changes in properties of the plasma and magnetic 
field, and are identified as statically advected tangential discontinuities (TDs), or propagating 
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rotational discontinuities (RDs). The TDs are characterized by small normal components of the 
magnetic field, large variations of magnetic field intensity and density jumps, while the RDs have 
large normal components of the magnetic field, but small variations of magnetic field intensity and 


density (Hudson 1970) 


By measurements of the magnetic field from Pioneer 6, Burlaga (1969) observed the disconti 


nuities in the magnetic-field direction with special emphasis on their distribution in time. Based 
on interplanetary field measurements made by the vector helium magnetometers onboard Pioneer 


10 and Pioneer 11, Tsurutani &; Smith (1979) investigated a possible dependence of the occurrence 


rate and the properties of the discontinuities on radial distance between 1 and 8.5 AU. Lepping V 
Behannon (1986) surveyed the data from the Mariner 10 primary mission to study the characteris¬ 


tics of the discontinuities in the interplanetary magnetic field at heliographic distances of 1.0, 0.72, 
and 0.46 AU, and found an j--- * 1 * ' 3 ^ 0 - 4 dependence for the daily average number of discontinuities per 
hour. With Ulysses magnetic field and plasma data obtained at radial distances ranging between 


1 and 5 AU from the Sun and at high heliographic latitudes, Tsurutani et al. (1994) discovered 


two regions where the occurrence rate of interplanetary discontinuities is high: in stream-stream 
interaction regions and in Alfven wave trains. To determine the normals of the discontinuities, 


Horbury et al. (2001) explored the discontinuities measured by three spacecraft WIND, IMP 8, and 


Geotail together with the solar wind velocity measured at Geotail, and obtained quite different dis¬ 


tributions of the discontinuity types. With magnetic field data from the ACE spacecraft, Vasquez 


et al. (2007) extended the survey of discontinuity properties to small spread angles of the held 


vectors across the discontinuity, and found that solar wind discontinuities are far more abundant at 


small than at large spread angles. By measurements from the WIND spacecraft, Wang et al. (2013) 


studied the intermittent structures in solar wind turbulence, which are identified as being mostly 
rotational discontinuities (RDs) and rarely tangential discontinuities (TDs) based on the technique 
described by Smith (1973) and |Tsurutani Sz Smith (1979). Paschmann et al. (2013) carried out 
a comprehensive study of directional discontinuties and Alfvenic fluctuations in the solar wind on 
the basis of Cluster data. 

So far, there are still debates regarding the origin and nature of the discontinuities in the 
solar wind. Since RDs appear as a compressed Alfven wave, nonlinear wave steepening has been 
suggested as the cause of its formation (Cohen &; Kulsrud|1974 Malar a Sz Elaouhr||1991 Tsurutani 


et al. 1994 Vasquez fe Hollweg|1996[ T~998a|b Tsurutani et al.|20~02 Marsch &; Verscharen|2011 ). By 
numerically calculating the evolution of an initially parallel-propagating, elliptically polarized wave 


train in a cold plasma, Cohen & Kulsrud (1974) were the first to investigate this possibility, and 


showed that this wave evolves to a constant-B solution with RDs that rotate the field by exactly 
180°. Vasquez &; Hollweg (1996) continued this study, but conducted a 1.5-D hybrid numerical 


simulation study of the evolution of obliquely propagating, linearly polarized Alfven wave trains. 
They found that large-amplitude cLB/Bq > 1 wave trains steepen and produce RDs which always 
rotate the held by < 180°. Vasquez & Hollweg (1998a) also presented 2.5-D numerical simulations 


of a small group of nonplanar Alfven waves to show the generation of imbedded RDs. It should 
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be noted that in their models the formation of RDs probably occurs relatively near the Sun where 
most Alfvenic fluctuations originate. 


There is another model suggesting that MHD turbulence dynamically generates these discon¬ 
tinuities as the solar wind flows outward. Recently, numerical simulations have been done to inves¬ 
tigate this assertion (Greco et al. 2008 2009, 2010; Servidio et al. 2011 Zhdankin et al. 2012a|b ). 
Through analyses of MHD simulation data, Greco et al. (2008) examined the relationship between 
discontinuities identified by classical methods, and coherent structures identified by using inter- 
mittency statistics. They found that the two methods produce remarkably similar distributions of 
waiting times, and in fact identify many of the same events. Greco et al. (2009) further examined 
the link between intermittent turbulence and MHD discontinuities, directly comparing simulations 
of MHD turbulence with statistical analysis from ACE solar wind data. Their results support the 
notion that some solar wind discontinuities are consequences of intermittent turbulence. In direct 


numerical simulations of MHD turbulence with an imposed uniform magnetic field, Zhdankin et al. 


(2012a) investigated the statistical properties of magnetic discontinuities, and concluded that the 
discontinuities observed in the solar wind can be reproduced by MHD turbulence. However, these 
works conducted statistical studies, and did not give a clear illustration of how discontinuities, 
especially RDs, are formed in MHD turbulence. 


In the present study, we utilize a compressible 3-D MHD model to illustrate and analyze 
the forming of RDs in the turbulence. By checking the magnetic field and plasma properties, it is 
found that RD is produced by the steepening of a moderate-amplitude Alfven wave with nonuniform 
propagating speed. The paper is organized as follows. In Section 2, a general description of the 
numerical MHD model is given. Section 3 describes the results of the numerical simulation and its 
analysis. Section 4 is reserved for the summary and discussion. 


2. NUMERICAL MHD MODEL 

The description of the plasma is given by compressible 3-D MHD, which involves a fluctu¬ 
ating flow velocity v(x,y, z,t), magnetic field h(x,y,z,t), density p(x,y, z,t), and temperature 
T(x,y, z,t). An uniform guide field Bq is assumed in the z—direction, so the total magnetic field 

is B = Bo + b. 

The MHD equations are written in the following non-dimensional form: 


dp 

dt 


+ V • pu = 0 


dpu 


+ V 


puu + (p+ -B 2 )I ~ BB 


= 0 


(1) 


dt 


( 2 ) 
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de 

dt 


+ V- 


u(e + p + ^B 2 ) - (u • B)B 


V • (B x n j), 


( 3 ) 


— + V • (uB - Bu) = ? ? V 2 B , (4) 

where 

e = \pv 2 + — + ^B 2 , j = V x B, (5) 

2 7—12 

which corresponds to the total energy density and current density, respectively. Here, p is the mass 
density; u = (v x , v y , v z ) are the x, y , and z—components of velocity; p is the thermal pressure; 
B denotes the magnetic field; t is time; 7 = 5/3 is the adiabatic index; and rj is the magnetic 
resistivity. 


Three independent parameters, an initial mean density po, a characteristic length L, and a 
characteristic plasma speed no = db/y/Anpo with 5b = ( 6 2 ) 1 / 2 , are used to normalize the MHD 
equations. Other variables are normalized by their combinations. The dimensionless numbers 
appearing in the equations are the Mach number M s = vq/c s , where c s = \fjpoJ~po is the sound 
speed, and the magnetic Reynolds number R m = vqL/ rj. Here, we take M s to be 0.5, consistent 
with the solar wind observations at 1 AU, and R m to be 1000, which is limited by the available 
spatial resolution. The uniform guide field Bq is two times of the fluctuating field | <5b |. 


We consider periodic boundary conditions in a cube with a side length of 2irL and a resolution 
defined by the number of grid points which is 512 3 , and run a simulation from an initial state with 
kinetic and magnetic energy per unit mass (v 2 ) = (b 2 ) = 1. The fluctuations initially populate an 
annulus in the Fourier A;—space such that 1 < k < 8 , with constant amplitude and random phases 
(Matthaeus et al. 1996 Dmitruk et al. 2004). The initial normalized cross helicity is set to be 0.9 
such that the primordial fluctuations are highly Alfvenic. The initial density and thermal pressure 
are set to be uniform. 


To solve the equations, we employ a splitting-based finite-volume numerical scheme. The fluid 
part is solved by the Godunov-type central scheme and the magnetic part by the constrained trans¬ 
port approach, in conjunction with the method called second-order Monotone Upstream Schemes 
for Conservation Laws (MUSCL) for reconstruction and with the approximate Riemann solvers of 
Harten-Lax-van Leer (HLL) for calculation of the numerical fluxes (Feng et al. 2011). The explicit 
second-order Runge-Kutta stepping with total variation diminishing is applied in time integration. 


3. NUMERICAL RESULTS 

Figure [T] shows the distribution of the current density in the z—direction J z in the x-z (Left) 
and x-y (Right) plane at t = 1.25. The arrows superposed on the images are the projections of 
the magnetic field vectors in the x-z (Left) and the x-y plane (Right). The black ellipses mark the 
region where the identified RD is formed. From this figure, we can see that a large-scale background 
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Fig. 1. Distribution of the z—component of the current density J z in an x-z (Left) and x-y 
(Right) plane at t = 1.25. Superposed by arrows are the projections of the magnetic field vectors 
in the x-z (Left) and the x-y plane (Right). The black ellipses mark the region where the identified 
RD is formed. 


magnetic field is clearly present in the direction. As a result of the well-known anisotropic 


behavior of magnetic field fluctuations in MHD with an imposed uniform guide field (Matthaeus 


et_al. 1996), current density structures preferentially align along the guide field direction, as shown 
in the left panel, and become much more varying in the perpendicular cross section, as shown in 
the right panel. Also, J z appears to be large in magnitude at the location where the identified RD 
is formed. 

In order to detect RD, we first seek the regions with large normalized partial variance of 
increments (PVI) of the magnetic field vector. PVI in 3-D space is defined as 


PVI (x,y,z) = 


B(x + Ax, y + Ay, z + A z) — B(x — Ax, y — Ay, z — Az) 


+ Ax, y + Ay, z + A z) — B(x — Ax, y — Ay, z — Az) | 2 } 


where Ax, Ay and A z is the grid increment in the x—, y— and z—direction, respectively. We 
first sample the magnetic field B, plasma velocity v, density p, and temperature T along a linear 
path through the region with PVI > 2, and then perform along that path the minimum variance 


analysis (MVA) (Sonnerup & Cahill 1967) using the magnetic field data to find the maximum 


variance direction ( L ), intermediate variance direction (M), minimum variance direction ( N ), and 
their corresponding eigenvalues Ai, A 2 , and A 3 . Finally, we check the 3-D plasma and field structure 
of the possible events. 


an 


Figure [2] shows the magnetic field and plasma parameters along the sampling interval for 
identified RD. In this figure, the magnetic field has been converted into the Alfven velocity Va, 
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Fig. 2.— Magnetic field and plasma parameters along the sampling interval for an identified RD. 
Vl, Vm-, and Vn denote the L,M,N components of the Alfven velocity (black lines) and plasma 
velocity (red lines) respectively as derived from minimum variance analysis (MVA). The Horizontal 
axis displays the coordinate s along the sampling path, with s = 0 at the RD point. 

by using Va = B/^/p . We can see that the sampling series of PVI is rather close to 0, except 
near the center. The large jumps of the Alfven velocity and plasma velocity mainly occur along 
the L direction, with Vl jumping from -1.0 to 1.0. The IV—components Vn of them are nearly 
constant, and are equal to 1.6, along with a nearly negligible jump of | Va |, the magnitude of 
Alfven speed. Also, the Alfven speed and plasma speed are nearly in agreement over the entire 
interval, including the jump across the RD itself. The positive correlation between them implies 
the RD propagation anti-parallel to the magnetic field. The density p, temperature T and total 
pressure Pt (which consists of the summed magnetic pressure and thermal pressure) all exhibit 
relatively constant traces throughout the interval. To be noted, the jump of Vl and the large value 
of Vn, together with the relatively slight change of Pt as well as p, corroborate that this event is a 
RD. 

Figure [3] exhibits the 3-D structure of the identified RD. The green lines in the left panel denote 
magnetic field lines and the yellow arrows in the right panel are plasma velocity vectors, which are 
converted into the de Hoffman-Teller (HT) frame. The light gray surface is the isosurface where 
PVI = 4, showing RD, and red, green as well as blue arrows display the x—, y—, and z—direction, 
respectively. This figure displays that in the HT frame, magnetic field and plasma velocity have 
evident normal components across the RD, and they both rotate by a certain angle. Also, the 
plasma velocity is field aligned on both sides of the discontinuity, in accordance with the Alfvenic 
nature of RD. The identified RD appears as a thin surface and its normal inclines to the z—axis. 

To see how the identified RD is formed, Figure [4] shows the isosurfaces of B x (Left) and V x 

















Fig. 3.— Three-dimensional structure of the identified RD. The green lines in the left panel 
denote magnetic field lines and the yellow arrows in the right panel are plasma velocity vectors, 
which are converted into de Hoffman-Teller (HT) frame. The light gray surface is the isosurface 
where PVI = 4, showing the RD, and red, green as well as blue arrows display the x—, y—, and 
z—direction, respectively. 



Fig. 4.— The isosurfaces of B x (Left) and V x (Right) at different moments in time, with the light 
gray surface being the isosurface where PVI = 3, and exhibiting the identified RD. 


















(Right) at different moments in time. These isosurfaces are associated with the Alfven wave as 
shown below. PVI = 3 is used here to exhibit the identified RD at these four moments, and is shown 
as the light gray surface. From this figure, it is notable that the mutual approaching of the two 
isosurfaces of B x ( V x ), which describes the steepening of the Alfven wave, leads to the formation 
of a RD. At t = 0.95, the isosurface with B x = —0.55 ( V x = —0.45) is relatively away from the 
isosurface with B x = 0.75 (V x = 0.45), and the transition of B x ( V x ) from B x = —0.55 ( V x = —0.45) 
to B x = 0.75 ( V x = 0.45) is gentle. The isosurface where PVI = 3 is small and thin. At t = 1.10, 
the isosurfaces with B x = —0.55 (14 = —0.45) and B x = 0.75 ( V x = 0.45) approach each other, 
and the transition between the two isosurfaces of B x ( V x ) becomes steep. Correspondingly, the 
isosurface where PVI = 3 grows. Afterwards, that is at t = 1.25, the two isosurfaces of B x (V x ) 
are close enough to each other, and the transition between the two isosurfaces of B x ( 14 ) becomes 
steeper. The isosurface where PVI = 3 grows larger and thicker, and the identified RD is fully 
grown. However, at t = 1.40, the isosurfaces with B x = —0.55 (14 = —0.45) is far away from the 
isosurface with B x = 0.75 (14 = 0.45), and the transition between the two isosurfaces of B x (V x ) 
becomes gentle again. As a result, the isosurface where PVI = 3 becomes small and thin. The 
identified RD starts to collapse. 

To understand the type of waves before the RD is formed and to investigate the process of 
the evolution of the isosurfaces mentioned above, Figure [5] exhibits the distribution of B x (the first 
row), 14 (the second row), and Va (the third row) in the neighborhood of the identified RD in the 
plane x-z at different moments in time. Superposed by arrows are the projections of the magnetic 
field (the first row), velocity field in the HT frame (the second row), and negative Alfven speed field 
(the third row) vectors in the x-z plane. The black ellipses mark the position where the identified 
RD is formed. Near the ellipses, the magnitudes of B x and 14 are almost same (they are set to 
have the same color scales so that the same color stands for the same value), and the directions of 
the in-plane projection of the B and V vectors are almost identical. Hence in this vicinity, we have 
V = B /y/~po (we recall that po ~ 1), which agrees with the polarity relations of an Alfven wave. In 
other words, the RD is detected in a neighbouring Alfvenic environment that apparently favours 
RD formation. 

Therefore, Va can be regarded as the propagation velocity of the structures relative to the 
location where the RD forms. Hence it is significant to trace the changes of the Alfven speed. In 
the third row in Figure [5] where the evolution of the Alfven speed is shown, we can see that at 
t = 0.95, there is a difference of Alfven speed across the black ellipse, which makes the layers with 
negative B x (blue in Figure [4]), propagate faster than that with positive B x (green in Figure [5]). At 
t = 1.10, there is the evidence for approaching and squeezing of these two layers. The difference 
of Alfven speed there remains. This will drive these two layers further to approach each other. 
At t = 1.25, their transition becomes sharp, and the difference of Alfven speed nearly disappears. 
This status continues until t = 1.40, when the transition becomes gentle as a result of the faster 
propagation of the layers with positive B x than that with negative B x . It is obvious that the 
difference of Alfven speed makes the Alfven wave steepen, a process which forms the identified RD. 
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Fig. 5.— Distribution of B x (the first row), V x (the second row), and Va (the third row) in a 
subzone of the x-z plane (which passes through the identified RD) at different moments in time. 
Superposed by arrows are the projections of the magnetic field (the first row), velocity field in the 
HT frame (the second row), and negative Alfven speed field (the third row) vectors in the x-z 
plane. The black ellipses mark the position where the identified RD is formed. 














































4. SUMMARY AND DISCUSSION 


In this study, we use a simulation of the decaying compressive 3-D MHD turbulence with an 
imposed uniform guide field as a test case to explore the formation of RDs in MHD turbulence. 
Motivated by solar wind observation at 1 AU, we consider a moderate fluctuation amplitude corre¬ 
sponding to 5b/B q = 0.5 and high Alfvenic correlations with the normalized cross helicity initially 
of 0.9. A case study is thus conducted to illustrate the origin of RDs in MHD turbulence. 

The numerical simulation shows the well-known anisotropic behavior of the turbulent MHD 
field, with the current density structures preferentially aligning along the guide field direction and 
scattering in the perpendicular plane. To detect the RDs in this simulated magnetofluid, we first 
seek the regions with large PVI, then conduct a MVA by sampling the parameters along a linear 
path, and finally check the 3-D structure of the possible RD events. One clear RD is identified with 
sharp field rotations and little variations of the magnetic field intensity as well as density. At the 
same time, in the HT frame, the plasma velocity is nearly in agreement with the Alfven speed, and 
is field aligned on both sides of the discontinuity, satisfying the Walen relation that expresses the 
Alfvenic nature of an RD. The normal direction of the identified RD inclines to the z—axis, and 
propagates anti-parallel to the guide field. 

The comprehensive information obtained by the simulation of the magnetic field and plasma 
parameters associated with the RD implies that the RD is produced by the steepening of the 
moderate-amplitude Alfven wave with nonuniform Alfven speed in the ambient turbulence. Before 
the RD is formed, the layers with negative B x smoothly transits to the layers with positive B x . 
However, there is a difference of Alfven speed across them, which makes the layers with negative B x 
chase after its counterpart with positive B x . As they are driven by the neigbouring turbulence to 
approach and squeeze each other, the transition between them undergoes further steepening until 
the difference of Alfven speed nearly disappears. At the same time, the identified RD is formed. 

In this work we investigated only one RD. Certainly, there are many RDs generated in the 
simulation of decaying MHD turbulence. It will be worthwhile to conduct statistical studies to see 
whether all RDs are produced by this nonlinear steepening of an Alfven wave. Meanwhile, TDs 
are also formed. It can be inspiring to investigate their generation mechanism in MHD turbulence. 
Furthermore, the parameters we adopted in the simulation (e.g. Mach number, cross helicity, plasma 
f3) may influence the forming of RDs or TDs. In the future, we plan to conduct a parameter study 
to investigate the possible effects induced by parameter variation on the generation of discontinuties 
in MHD turbulence. 


This work is supported by NSFC under contract Nos. 41304133, 41474147, 41231069, 41222032, 
41174148, 41421003, and 41204105. The numerical calculation has been completed on computing 
system at Peking University. 
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